﻿/* A C-program for MT19937, with initialization improved 2002/1/26.
   Coded by Takuji Nishimura and Makoto Matsumoto.

   Before using, initialize the state by using init_genrand(seed)  
   or init_by_array(init_key, key_length).

   Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,
   All rights reserved.                          

   Redistribution and use in source and binary forms, with or without
   modification, are permitted provided that the following conditions
   are met:

     1. Redistributions of source code must retain the above copyright
        notice, this list of conditions and the following disclaimer.

     2. Redistributions in binary form must reproduce the above copyright
        notice, this list of conditions and the following disclaimer in the
        documentation and/or other materials provided with the distribution.

     3. The names of its contributors may not be used to endorse or promote 
        products derived from this software without specific prior written 
        permission.

   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
   "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
   LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
   A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT OWNER OR
   CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
   EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
   PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
   PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
   LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
   NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
   SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 * */

/* 这是 32位的 梅森旋转算法 （Mersenne Twister pseudorandom number generator）
 * 参考了 Takuji Nishimura 和 Makoto Matsumoto 的 MT19937ar.sep (2001/1/26 version) 实现
 * 对具体实现进行了一些修改
 * 
 * 梅森旋转算法 主页：
 * http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/MT2002/emt19937ar.html
 * */

using System;

namespace HS.Math
{
    class MTRandom
    {
        private const int N = 624;
        private const int M = 397;
        private const uint MATRIX_A = 0x9908b0dfU; //魔数

        private const uint MATRIX_INIT = 0x6C078965u; //1812433253U
        private const uint MATRIX_FORE = 0x19660Du; //1664525U
        private const uint MARTIX_BACK = 0x5d588b65u; //1566083941U

        private const uint OUT_SECOND = 0x9d2c5680u;
        private const uint OUT_THIRD = 0xefc60000u;

        private const uint UPPER_MASK = 0x80000000U; //最重要的 1 个 bit
        private const uint LOWER_MASK = 0x7fffffffU; //最不重要的 31 个 bit

        private const uint MBS = 0x80000000u; //0x01u << 31;
        private const uint MAX = 0xFFFFFFFFu;

        private uint[] mt = new uint[N]; //输出的数据表
        private uint index = N + 1; //无论采用何种方式创建类对象，创建后第一次输出数据前都会混洗一次

        private void Init(uint seed)
        {
            //初始化 mt
            mt[0] = seed;
            for (index = 1; index < N; index++)
            {
                mt[index] = (MATRIX_INIT * (mt[index - 1] ^ (mt[index - 1] >> 30)) + index);
            }
        }

        private void Init(uint[] seed)
        {
            int i, j, k;
            i = 1; //index of mt
            j = 0; //index of init
            k = (N > seed.Length ? N : seed.Length);
            for (; k > 0; k--)
            {
                mt[i] = (mt[i] ^ ((mt[i - 1] ^ (mt[i - 1] >> 30)) * MATRIX_FORE)) + seed[j] + (uint)j;  /* non linear */
                i++; j++;

                //回到初始位置，防止下标溢出
                if (i >= N) { mt[0] = mt[N - 1]; i = 1; }
                if (j >= seed.Length) j = 0;
            }
            for (k = N - 1; k > 0; k--)
            {
                mt[i] = (mt[i] ^ ((mt[i - 1] ^ (mt[i - 1] >> 30)) * MARTIX_BACK)) - (uint)i; /* non linear */
                i++;

                //回到初始位置，防止下标溢出
                if (i >= N) { mt[0] = mt[N - 1]; i = 1; }
            }

            mt[0] = MBS;
        }

        public MTRandom()
        {
            //入口掩码
            uint next;
            next = (uint)DateTime.Now.ToBinary();
            next += (uint)DateTime.Now.Millisecond;

            //生成一组初始 mt[]
            Init(next);
        }

        public MTRandom(uint seed)
        {
            //入口掩码
            uint next = seed;
            next ^= (uint)DateTime.Now.ToBinary();
            next += (uint)DateTime.Now.Millisecond;

            //生成一组初始 mt[]
            Init(next);
        }

        public MTRandom(uint[] seed)
        {
            //int i, j, k;

            //入口掩码
            uint next;
            next = (uint)DateTime.Now.ToBinary();
            next += (uint)DateTime.Now.Millisecond;

            //生成一组初始 mt[]
            Init(next);
            Init(seed);
        }

        //返回值范围：[0, MAX]
        internal uint Next()
        {
            int i;
            uint next;
            uint[] mag01 = new uint[2] { 0x0u, MATRIX_A };

            //数组的数据第一次使用前也执行一次
            //如果数组的数据都用了一遍，也重新生成数据
            //N 个数据为一个周期 
            if (index >= N)
            {
                //第 [0, N - M) 个数据 使用第一种方式处理
                for (i = 0; i < N - M; i++)
                {
                    next = (mt[i] & UPPER_MASK) | (mt[i + 1] & LOWER_MASK);
                    mt[i] = mt[i + M] ^ (next >> 1) ^ mag01[next & 0x1u];
                }

                //第 [N - M, N - 1) 个数据 使用第二种方式处理
                for (; i < N - 1; i++)
                {
                    next = (mt[i] & UPPER_MASK) | (mt[i + 1] & LOWER_MASK);
                    mt[i] = mt[i + M - N] ^ (next >> 1) ^ mag01[next & 0x1u];
                }

                //第 N - 1 个数据 使用第三种方式处理
                next = (mt[N - 1] & UPPER_MASK) | (mt[0] & LOWER_MASK);
                mt[N - 1] = mt[M - 1] ^ (next >> 1) ^ mag01[next & 0x1u];

                //重置索引
                index = 0;
            }

            //获取数组中轮到输出的数据
            next = mt[index++];

            //混乱数据
            next ^= (next >> 11);
            next ^= (next << 7) & OUT_SECOND;
            next ^= (next << 15) & OUT_THIRD;
            next ^= (next >> 18);

            return next;
        }

        internal uint NextUInt(uint minValue, uint maxValue)
        {
            uint next;
            double temp;

            next = Next();

            //按比例转换成对应范围
            temp = next / (MAX + 1.0);
            temp = temp * (maxValue - minValue) + minValue;

            return (uint)temp;
        }

        //与 .NET 不同，支持输出负数
        internal int NextInt(int minValue, int maxValue)
        {
            uint next;
            double temp;

            next = Next();

            //按比例转换成对应范围
            temp = next / (MAX + 1.0);
            temp = temp * (maxValue - minValue) + minValue;

            return (int)temp;
        }

        //返回值范围：[0, 1)
        internal double NextDouble()
        {
            uint next;
            double temp;

            next = Next();

            //按比例转换成对应范围
            temp = next / (MAX + 1.0);

            return temp;
        }
    }
}